{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine Learning Exercise 3 - Multi-Class Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook covers a Python-based solution for the third programming exercise of the machine learning class on Coursera.  Please refer to the [exercise text](https://github.com/jdwittenauer/ipython-notebooks/blob/master/exercises/ML/ex3.pdf) for detailed descriptions and equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "For this exercise we'll use logistic regression to recognize hand-written digits (0 to 9).  We'll be extending the implementation of logistic regression we wrote in exercise 2 and apply it to one-vs-all classification.  Let's get started by loading the data set.  It's in MATLAB's native format, so to load it in Python we need to use a SciPy utility."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'X': array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "        [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "        [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "        ..., \n",
       "        [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "        [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "        [ 0.,  0.,  0., ...,  0.,  0.,  0.]]),\n",
       " '__globals__': [],\n",
       " '__header__': 'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011',\n",
       " '__version__': '1.0',\n",
       " 'y': array([[10],\n",
       "        [10],\n",
       "        [10],\n",
       "        ..., \n",
       "        [ 9],\n",
       "        [ 9],\n",
       "        [ 9]], dtype=uint8)}"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat\n",
    "%matplotlib inline\n",
    "\n",
    "data = loadmat('data/ex3data1.mat')\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000L, 400L), (5000L, 1L))"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data['X'].shape, data['y'].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, we've got our data loaded.  The images are represented in martix X as a 400-dimensional vector (of which there are 5,000 of them).  The 400 \"features\" are grayscale intensities of each pixel in the original 20 x 20 image.  The class labels are in the vector y as a numeric class representing the digit that's in the image.\n",
    "\n",
    "The exercise code in MATLAB has a function provided to visualize the hand-written digits.  I'm not going to reproduce that in Python, but there's an illustration in the exercise PDF if one is interested in seeing what the images look like.  We're going to move on to our logistic regression implementation.\n",
    "\n",
    "The first task is to modify our logistic regression implementation to be completely vectorized (i.e. no \"for\" loops).  This is because vectorized code, in addition to being short and concise, is able to take advantage of linear algebra optimizations and is typically much faster than iterative code.  However if you look at our cost function implementation from exercise 2, it's already vectorized!  So we can re-use the same implementation here.  Note we're skipping straight to the final, regularized version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    return 1 / (1 + np.exp(-z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cost(theta, X, y, learningRate):\n",
    "    theta = np.matrix(theta)\n",
    "    X = np.matrix(X)\n",
    "    y = np.matrix(y)\n",
    "    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))\n",
    "    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))\n",
    "    reg = (learningRate / 2 * len(X)) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))\n",
    "    return np.sum(first - second) / (len(X)) + reg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we need the function that computes the gradient.  Again, we already defined this in the previous exercise, only in this case we do have a \"for\" loop in the update step that we need to get rid of.  Here's the original code for reference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gradient_with_loop(theta, X, y, learningRate):\n",
    "    theta = np.matrix(theta)\n",
    "    X = np.matrix(X)\n",
    "    y = np.matrix(y)\n",
    "    \n",
    "    parameters = int(theta.ravel().shape[1])\n",
    "    grad = np.zeros(parameters)\n",
    "    \n",
    "    error = sigmoid(X * theta.T) - y\n",
    "    \n",
    "    for i in range(parameters):\n",
    "        term = np.multiply(error, X[:,i])\n",
    "        \n",
    "        if (i == 0):\n",
    "            grad[i] = np.sum(term) / len(X)\n",
    "        else:\n",
    "            grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])\n",
    "    \n",
    "    return grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our new version we're going to pull out the \"for\" loop and compute the gradient for each parameter at once using linear algebra (except for the intercept parameter, which is not regularized so it's computed separately).  To follow the math behind the transformation, refer to the exercise 3 text.\n",
    "\n",
    "Also note that we're converting the data structures to NumPy matrices (which I've used for the most part throughout these exercises).  This is done in an attempt to make the code look more similar to Octave than it would using arrays because matrices automatically follow matrix operation rules vs. element-wise operations, which is the default for arrays.  There is some debate in the community over wether or not the matrix class should be used at all, but it's there so we're using it in these examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gradient(theta, X, y, learningRate):\n",
    "    theta = np.matrix(theta)\n",
    "    X = np.matrix(X)\n",
    "    y = np.matrix(y)\n",
    "    \n",
    "    parameters = int(theta.ravel().shape[1])\n",
    "    error = sigmoid(X * theta.T) - y\n",
    "    \n",
    "    grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)\n",
    "    \n",
    "    # intercept gradient is not regularized\n",
    "    grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)\n",
    "    \n",
    "    return np.array(grad).ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've defined our cost and gradient functions, it's time to build a classifier.  For this task we've got 10 possible classes, and since logistic regression is only able to distiguish between 2 classes at a time, we need a strategy to deal with the multi-class scenario.  In this exercise we're tasked with implementing a one-vs-all classification approach, where a label with k different classes results in k classifiers, each one deciding between \"class i\" and \"not class i\" (i.e. any class other than i).  We're going to wrap the classifier training up in one function that computes the final weights for each of the 10 classifiers and returns the weights as a k X (n + 1) array, where n is the number of parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "def one_vs_all(X, y, num_labels, learning_rate):\n",
    "    rows = X.shape[0]\n",
    "    params = X.shape[1]\n",
    "    \n",
    "    # k X (n + 1) array for the parameters of each of the k classifiers\n",
    "    all_theta = np.zeros((num_labels, params + 1))\n",
    "    \n",
    "    # insert a column of ones at the beginning for the intercept term\n",
    "    X = np.insert(X, 0, values=np.ones(rows), axis=1)\n",
    "    \n",
    "    # labels are 1-indexed instead of 0-indexed\n",
    "    for i in range(1, num_labels + 1):\n",
    "        theta = np.zeros(params + 1)\n",
    "        y_i = np.array([1 if label == i else 0 for label in y])\n",
    "        y_i = np.reshape(y_i, (rows, 1))\n",
    "        \n",
    "        # minimize the objective function\n",
    "        fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)\n",
    "        all_theta[i-1,:] = fmin.x\n",
    "    \n",
    "    return all_theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few things to note here...first, we're adding an extra parameter to theta (along with a column of ones to the training data) to account for the intercept term.  Second, we're transforming y from a class label to a binary value for each classifier (either is class i or is not class i).  Finally, we're using SciPy's newer optimization API to minimize the cost function for each classifier.  The API takes an objective function, an initial set of parameters, an optimization method, and a jacobian (gradient) function if specified.  The parameters found by the optimization routine are then assigned to the parameter array.\n",
    "\n",
    "One of the more challenging parts of implementing vectorized code is getting all of the matrix interactions written correctly, so I find it useful to do some sanity checks by looking at the shapes of the arrays/matrices I'm working with and convincing myself that they're sensible.  Let's look at some of the data structures used in the above function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000L, 401L), (5000L, 1L), (401L,), (10L, 401L))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rows = data['X'].shape[0]\n",
    "params = data['X'].shape[1]\n",
    "\n",
    "all_theta = np.zeros((10, params + 1))\n",
    "\n",
    "X = np.insert(data['X'], 0, values=np.ones(rows), axis=1)\n",
    "\n",
    "theta = np.zeros(params + 1)\n",
    "\n",
    "y_0 = np.array([1 if label == 0 else 0 for label in data['y']])\n",
    "y_0 = np.reshape(y_0, (rows, 1))\n",
    "\n",
    "X.shape, y_0.shape, theta.shape, all_theta.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These all appear to make sense.  Note that theta is a one-dimensional array, so when it gets converted to a matrix in the code that computes the gradient, it turns into a (1 X 401) matrix.  Let's also check the class labels in y to make sure they look like what we're expecting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10], dtype=uint8)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.unique(data['y'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make sure that our training function actually runs, and we get some sensible outputs, before going any further."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -5.79312170e+00,   0.00000000e+00,   0.00000000e+00, ...,\n",
       "          1.22140973e-02,   2.88611969e-07,   0.00000000e+00],\n",
       "       [ -4.91685285e+00,   0.00000000e+00,   0.00000000e+00, ...,\n",
       "          2.40449128e-01,  -1.08488270e-02,   0.00000000e+00],\n",
       "       [ -8.56840371e+00,   0.00000000e+00,   0.00000000e+00, ...,\n",
       "         -2.59241796e-04,  -1.12756844e-06,   0.00000000e+00],\n",
       "       ..., \n",
       "       [ -1.32641613e+01,   0.00000000e+00,   0.00000000e+00, ...,\n",
       "         -5.63659404e+00,   6.50939114e-01,   0.00000000e+00],\n",
       "       [ -8.55392716e+00,   0.00000000e+00,   0.00000000e+00, ...,\n",
       "         -2.01206880e-01,   9.61930149e-03,   0.00000000e+00],\n",
       "       [ -1.29807876e+01,   0.00000000e+00,   0.00000000e+00, ...,\n",
       "          2.60651472e-04,   4.22693052e-05,   0.00000000e+00]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_theta = one_vs_all(data['X'], data['y'], 10, 1)\n",
    "all_theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're now ready for the final step - using the trained classifiers to predict a label for each image.  For this step we're going to compute the class probability for each class, for each training instance (using vectorized code of course!) and assign the output class label as the class with the highest probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def predict_all(X, all_theta):\n",
    "    rows = X.shape[0]\n",
    "    params = X.shape[1]\n",
    "    num_labels = all_theta.shape[0]\n",
    "    \n",
    "    # same as before, insert ones to match the shape\n",
    "    X = np.insert(X, 0, values=np.ones(rows), axis=1)\n",
    "    \n",
    "    # convert to matrices\n",
    "    X = np.matrix(X)\n",
    "    all_theta = np.matrix(all_theta)\n",
    "    \n",
    "    # compute the class probability for each class on each training instance\n",
    "    h = sigmoid(X * all_theta.T)\n",
    "    \n",
    "    # create array of the index with the maximum probability\n",
    "    h_argmax = np.argmax(h, axis=1)\n",
    "    \n",
    "    # because our array was zero-indexed we need to add one for the true label prediction\n",
    "    h_argmax = h_argmax + 1\n",
    "    \n",
    "    return h_argmax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the predict_all function to generate class predictions for each instance and see how well our classifier works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy = 97.58%\n"
     ]
    }
   ],
   "source": [
    "y_pred = predict_all(data['X'], all_theta)\n",
    "correct = [1 if a == b else 0 for (a, b) in zip(y_pred, data['y'])]\n",
    "accuracy = (sum(map(int, correct)) / float(len(correct)))\n",
    "print 'accuracy = {0}%'.format(accuracy * 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Almost 98% isn't too bad!  That's all for exercise 3.  In the next exercise, we'll look at how to implement a feed-forward neural network from scratch."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
